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MIGRATION AND GROWTH OF GIANT PLANETS IN SELF-GRAVITATING DISKS WITH VARIED 

THERMODYNAMICS 
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ABSTRACT 

We report on the results of novel global high-resolution three-dimensional simulations of disk-planet in- 
teraction which incorporate simultaneously realistic radiation physics and the self-gravity of the gas, as well 
as allowing the planet to move. We find that thermodynamics and radiative physics have a remarkable ef- 
fect on both migration and accretion of Jupiter mass planets. In simulations with radiative transfer adopting 
flux-limited diffusion, inward migration can be decreased by about 30% relative to the isothermal case, while 
in adiabatic runs migration nearly shuts off after a few tens of orbits. Migration varies because the relative 
strength of the inner and outer spiral perturbations is affected by thermodynamics, thus changing the net torque 
acting on the planet. Mass accretion rates on the planet can be reduced by more than an order of magnitude 
going from isothermal to radiative transfer and adiabatic simulations. A circumplanetary disk always forms 
except in adiabatic runs. With radiative transfer the disk is sub-keplerian (y m thkep ~ 0.7) owing to significant 
pressure support. We discuss the effect of circumplanetary disk structure on the drift of embedded dust grains 
and planetesimals and thus on the formation of the rocky satellites of giant planets. 

Subject headings: Hydrodynamics - planetary systems: formation - planetary systems: protoplanetary disks 



1. INTRODUCTION 

In the conventional picture, giant planets form via a two- 
stage process, known as core-accretion. In such a model, first 
a massive rocky core is assembled via gravitational accumu- 
lation of planetesimals in the protoplanetary disk and then the 
core begins to accrete the surrounding gas once it has grown 
massive enough (Pollack et al. 1996; Ida & Lin 2004; Alib- 
ert et al. 2005). The growing planet excites density waves as 
it moves through the disk. Such waves exert a torque on the 
planet (Lin & Papaloizou 1993) whose net effect is to extract 
angular momentum from its orbit for standard protosolar neb- 
ula models (Ward 1997). 

Different regimes of migration have been identified by nu- 
merous studies depending on the planet's mass and on the im- 
portance of co-rotation torques (see Papaloizou et al. 2007, for 
a review). In particular, when the planet has a mass compara- 
ble to Saturn or larger the interaction with the disk becomes 
markedly nonlinear (type II migration) and the planet carves 
a gap in the disk as a result of the planetary tide (Lin & Pa- 
paloizou 1993; Crida et al. 2007). Then, the planet decouples 
dynamically from the disk and migrates inward at a pace de- 
termined by the local viscous timescale (Lin & Papaloizou 
1986). 

A major assumption of almost all simulations published so 
far is that the disk is locally isothermal. This means that any 
heating is immediately radiated away, which would only be 
true if the disk was optically thin. The planet moves super- 
sonically in the disk, shock-heating the gas, and gas accretion 
onto the planet generates compressional heating. With typical 
densities 10 _1 'g/cm 3 , the midplane of a minimum solar nebula 
disk is indeed optically thick (r w 10) (d' Angelo et al. 2003). 
Recently, Paardekooper & Mellema (2006, 2008, hereafter 
PM06; PM08) have used both 2D and 3D adaptive mesh re- 
finement simulations with radiative transfer modeled via flux- 
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limited diffusion finding that for low mass planets both inward 
migration and gas accretion can be strongly suppressed. Klahr 
& Kley (2006, hereafter KK06) studied Jupiter-sized planets 
with a static 3D grid code and flux-limited diffusion. They 
did not report significant differences on migration compared 
with the isothermal case but found the structural evolution of 
the circumplanetary gas distribution to be strongly affected 
by radiation physics as noticed earlier in the 2D nested-grid 
calculations by d' Angelo et al. (2003). 

Yet, even these recent simulations lack several important in- 
gredients. First, they do not treat self-consistently the dynam- 
ics of the disk, planet and star; the planet and the star cannot 
move, and in some cases a gap is introduced already at the 
beginning of the simulation (KK06). Second, except for the 
2D simulations of d' Angelo et al. (2003), they adopt inviscid 
disks. Recently, Edgar (2007) found that giant planet migra- 
tion in an isothermal viscous disk does not obey the standard 
type II regime. A deep gap is never produced and hence mi- 
gration does not proceed on the viscous timescale. Finally, 
all these simulations neglect the self-gravity of the gas, which 
is required when simulating freely moving planets and might 
affect disk torques (Baruteau & Masset 2008). Self-gravity 
has been previously included only in a few isothermal calcu- 
lations (Nelson & Benz 2003a,b; Lufkin et al. 2004). 

In this Letter we present the first high-resolution three- 
dimensional hydrodynamical simulations of the interaction 
between a massive, Jupiter-sized planet and a surrounding vis- 
cous protoplanetary disk that include simultaneously radiative 
transfer, shock heating, self-gravity of the gas and fully self- 
consistent dynamics. We study migration, mass flow towards 
the planet and the circumplanetary gas distribution, exploring 
also the implications on the formation of satellites of giant 
planets. We compare the results with those obtained for lo- 
cally isothermal disks as well as other simulations with sim- 
plified disk thermodynamics. 

2. INITIAL CONDITIONS AND SIMULATIONS 

We have run different simulations with varying disk masses, 
resolution and thermodynamics (see Table 1). Our simu- 
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Simulations parameters 
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k The Table shows the accretion and migration rates for the various runs. The 
legend for the run names is as follows: IsoT refers to Isothermal runs, FLD 
to runs with the Flux Limited Diffusion approximation, NIsoT are adiabatic 
runs with 7 = 6/5 and Adia are adiabatic runs with 7 = 7/5. 100K corresponds 
to 100 000 particles, 200K to 200 000 particles and 1M to 1 Million particles. 
The runs are normally done with a planet gravitationnal softening of R/y/5, 
where R« =0.35 AU is the Hill radius for a Jovian planet, except for those 
marked with * for which the softening is equal to R# . Accretion rates are 
given in Mj yr and migration rates in AU yr~' . Mean rates are computed 
by averaging over time over the full extent of the simulation. Partial rates 
are computed after a number of orbits equivalent to the maximum number 
of orbits for the FLD runs, namely after 20.6 orbits for model Ml and 27.9 
orbits for model M2. When the accretion rate is 0, it means that the mass 
inside the Hill radius is below the SPH mass resolution. 
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FIG. 1. — Temperature map with overplotted midplane density contours of 
run FLD1M2 after 2.9 orbits. The shock along the spiral arms triggered by 
the planet is evident. The planet is seen as a hot spot in the disk. 

lations include the self-gravity of the gas (Nelson & Benz 
2003a,b; Lufkin et al. 2004). We place a 1 Jupiter mass 
(IMj) planet on a circular orbit at 5 AU for model Ml (resp. 
5.4 AU for model M2) from a 1M© mass star in an initially 
axisymmetric disk. The planet and the star are represented 
by softened N-Body particles free to move in response to 
their mutual gravity and that of the disk. There is no gap 
initially. Mass is allowed to accumulate around the planet 
since we do not use a prescription for accretion. Our ini- 
tial disk extends from 2 to 20 AU in model Ml and from 
1 to 25 AU in model M2. The surface density profile is 
£ = S(5 AU)(r/5 AU)~ LS . We adopt an initial surface density 
of S(5 AU) = 75 g cm" 2 (resp. 150 g cm" 2 ), which yield a total 
disk mass of M d i s k = 0.004 M Q for model Ml and 0.01 M Q for 
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FIG. 2. — Semi-major axis of the planet versus time (left) and azimufhally 
averaged radial torques exerted by the disk on the planet (right) as a function 
of the distance from the star for runs IsoT100K2, FLD100K2, NIsoT100K2 
and Adial00K2 after 20 orbits. Note that, in the adiabatic case, the inner and 
outer torques nearly balance. The legend for the different curves is in the 
figure. The position of the planet is shown with symbols of different colors 
depending on the run (right panel). 

model M2. The vertical density profile is initially gaussian. 
We assume the same disk temperature profile as in KK06 for 
the sake of comparison, T = T^r' 1 (T is independent on z). 
The pressure scale-height is initially chosen to be H/r = 0.05 
and implies Tq= 102 K. 

As for disk thermodynamics, we consider the following dif- 
ferent cases: (1) locally isothermal (IsoT) runs, in which the 
temperature of individual particles is kept constant over time; 
(2) adiabatic (Adia) and nearly isothermal (NIsoT) runs, in 
which the gas is evolved adiabatically assuming an adiabatic 
index of, respectively, 7 = 7/5 and 7 = 6/5 and solving the in- 
ternal energy equation including shock heating; (3) radiative 
transfer (FLD) runs, in which we solve the energy equation 
with added flux-limited diffusion and blackbody radiation at 
the disk edge to model the thermal energy flow in the disk, as 
described in Mayer et al. (2007). We adopt realistic opacities 
by d'Alessio et al. (1997). 

The simulations are performed with the parallel 3D 
Tree+SPH code GASOLINE (Wadsley et al. 2004), which 
includes adaptive multiple timesteps. In radiative transfer 
runs radiative timesteps can become very short so that only a 
few tens of orbits could be explored (see also PM06; PM08). 
We use the standard Monaghan viscosity with coefficients 
a = 1 and (3 = 2 and with the Balsara switch to reduce 
viscosity significantly in shear flows (Balsara 1995). We 
find an average equivalent effective Shakura-Sunyaev ass 
parameter ~ 10" 2 , although this varies appreciably with time 
and location in the disk. For comparison, protoplanetary 
disks are expected to have ass between 10" 4 and 10" 1 as a 
result of MHD turbulence (Papaloizou & Nelson 2003). 

3. RESULTS 

We find that both the orbital evolution and the properties 
of the gas flow around the planet depend significantly on disk 
thermodynamic s . 

3.1. Migration 

Our reference case for the study of migration is repre- 
sented by the locally isothermal runs. These have been widely 
used in the literature, although only rarely in the past self- 
gravity was included and/or the planet was allowed to move 
(Papaloizou et al. 2007; Bate et al. 2003; Nelson & Benz 
2003a,b). The planet carves a gap quite slowly compared to 
simulations of inviscid disks (e.g. Nelson et al. 2003a). How- 
ever, the drop in the surface density near the planet reaches 
more than order order of magnitude after 190 orbits, in agree- 
ment with the results of other three-dimensional grid-based 
and SPH simulations of Jupiter-mass planets in disk with 
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moderate viscosity (de Val-Borro et al. 2006). At this late 
stage the migration timescale is ~ 3 x 10 4 years, somewhat 
faster than the timescale of ~ 8 x 10 4 years typically inferred 
for Jupiter mass planets evolving in inviscid disks on a fixed 
orbit (e.g. Papaloizou et al. 2007; KK06; Edgar 2007). A 
speed-up of the migration rate by a factor 2-3 for disks with 
a viscosities ~ 10 has already been reported in the literature 
(d'Angelo et al. 2003; Edgar 2007). For the first few tens of 
orbits, when the gap is barely opened, the planet excites strong 
spiral modes since it is still well coupled with the disk and the 
migration rate is much faster than in the late stage (see par- 
tial migration rate in Table 1). Self-gravity also should lead 
to an acceleration of migration before the planet has opened 
a gap (Baruteau & Masset 2008). The mean migration rate 
increases with increasing disk mass (Table 1), in agreement 
with Edgar (2007), who also studied viscous isothermal disks 
that do not open deep gaps. 

In flux-limited diffusion (FLD) runs the planet has a signif- 
icant effect on the temperature structure of the disk around it. 
The triggered spiral arms are associated with a shock since 
the planet is moving supersonically. The temperature along 
the spiral arms thus rises significantly because the disk mid- 
plane is optically thick (Figure 1) an effect which is obviously 
absent in isothermal runs. Flux-limited diffusion runs are fol- 
lowed up to a point when they are still in the weak gap regime, 
i.e. for no more than 30 orbits. The gap is somewhat weaker 
in these runs because the disk is heated by the spiral shocks 
in the region around the planet, becoming thicker and thus 
making it harder to satisfy the gap opening condition. 

Despite the shallower gap, in FLD runs migration is slower 
than in the isothermal case, by 30% for the case of the low 
disk mass, and by nearly 15% for the case of the low mass 
disk (Figure 2, left panel and Table 1). The slower migration 
is due to an effect already seen and explained in PM08 for 
deeply embedded, low-mass planets. It is due to the change 
of the balance between the outer and inner torque acting on 
the planet as the thermodynamics is varied. In particular, we 
find that, the net torque is more negative in the isothermal runs 
relative to the FLD ones (Figure 2, right panel). The variation 
of torque balance is due to a variation in the structure of the 
spiral perturbations excited by the planet. 

Torques are even more affected by the different tempera- 
ture and density evolution in the adiabatic runs (Figure 2, left 
panel). In these runs the gas along the spiral shocks rises its 
temperature by up to a factor of 5 and adiabatically expands 
away from the midplane. As a result, after 50 orbits the gas 
surface density in the disk midplane around the planet has de- 
creased by a factor 12.5 in the run with 7 = 7/5 and by 8.3 
in the run with 7 = 6/5. Because of the decreased gas surface 
density, the spiral arms around the planet weaken consider- 
ably with time, and so do the torques exerted by them. Then, 
after a few tens of orbits, migration nearly shuts off (Figure 2, 
right panel). The planet is never able to open a gap because 
the radius of the planet's Roche lobe is smaller than the local 
disk semi-thickness. 

3.2. Gas accretion and planetary growth 

Since we do not include explicitly accretion onto the planet 
in our simulations we calculate an hypothetical accretion rate 
by measuring the mass of gas that accumulates within the Hill 
radius of the planet (Rh = 0.35 AU) as a function of time. The 
accretion rate thus defined is found to be strongly dependent 
on disk thermodynamics as well as on the disk mass, and, not 
surprisingly, increases moderately as the gravitational soften- 




FlG. 3. — Azimuthally averaged temperature maps of circumplanetary gas 
with superimposed density countours forFLDlM2 (left) and IsoTlM2 (right) 
runs after 2.9 orbits. The colorbar gives the temperature on a logarithmic 
scale in K. 

ing of the planet is decreased (see Table 1). 

In isothermal Ml runs the accretion rate averaged over 
nearly 200 orbits is 5 x 10~ 5 My yr" 1 . in agreement with pre- 
vious estimates obtained for similar disk masses and effective 
a viscosities (e.g. Bate et al. 2003; d'Angelo et al. 2003). 
KK06 compute mass accretion rates about 3 times smaller for 
planets in an isothermal disk but their disk is inviscid and they 
start with a gap, both differences going in the direction of low- 
ering somewhat the mass accretion. 

As seen in PM08 for the case of low mass planets, planetary 
accretion in non-isothermal disks can be reduced compared 
to isothermal disks. This is because the gas is compression- 
ally heated as it flows towards the planet, and the resulting 
pressure gradient opposes further gas inflow. The temper- 
ature structure around the planet in one of the FLD runs is 
shown in Figure 3. For our lowest disk mass (Ml) the accre- 
tion rate over the first 20 orbits is ~ 2 x 10~ 5 My yr" 1 in the 
FLD run, 20 times lower than for the isothermal during the 
same part of the orbital evolution (Table 1), and comparable 
with what KK06 find for their radiative transfer runs starting 
with or without a gap. Similar rates were also obtained by 
d'Angelo et al. (2003) for their radiative viscous models. The 
difference between isothermal and FLD runs is smaller in the 
case of the largest disk mass (M2). In the FLD runs with in the 
M2 disk, which have the smallest softening, the planet should 
be able to accrete a few Jupiter masses over ~ 10 4 years, a 
time comparable to the migration rate. 

In adiabatic disks accretion, like migration, essentially 
shuts off after a few tens of orbits (Table 1). Due to the 
absence of cooling compressional heating in this case is in- 
deed more effective. Moreover, heating due to spiral shocks 
raises the temperature of the gas even at some distance from 
the planet, which further stifles accretion. 

3.3. The circumplanetary disk 

A circumplanetary disk begins to form already during the 
first few orbits in both the isothermal and the radiative trans- 
fer runs (Figure 3). The disk is thicker in the radiative trans- 
fer runs due to higher thermal pressure. It is sub-keplerian 
in both cases. Yet, there is a fairly extended region, between 
0.15 and 0.35 AU from the planet, within which the gas az- 
imuthal velocity is close to 70% of the local keplerian velocity 
(Figure 4, left panel). In the adiabatic runs the gas becomes 
hot enough to form a diffuse pressure supported atmosphere 
around the planet rather than a dense circumplanetary disk. In 
the latter case the gas is markedly sub-keplerian (Figure 4, 
left panel). In reality this is an extreme configuration that 
could occur only if cooling is completely suppressed by some 
stronger heating mechanism. A protoplanetary disk strongly 
irradiated by neighboring OB stars in a dense cluster might 
present such extreme conditions at large radii, R > 10 AU. 
(Alexander et al. 2006). 
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FIG. 4. — Gas azimuthal velocity divided by the keplerian velocity (left) and 
grain size for which T s /T m -t=\l(2 n), which correpond to the fastest possible 
migration towards pressure maxima (right) for runs IsoT100K2, FLD100K2, 
NIsoT100K2 and Adial00K2 after 27.9 orbits. The legend for the different 
curves is in the right panel. 

The velocity, density and temperature of the gas flow in 
the circumplanetary disk can determine whether rocky satel- 
lites of giant planets may form or not. Based on the strong 
sub-keplerianity of the gas around their Jupiter-mass planet, 
KK06 concluded that satellites could not form in their radia- 
tive transfer simulations. 

In Figure 4, right panel, we show the particle size for which 
T s /T or b=l/(2 7r) where T s is the stopping time in the Epstein 
regime and T ol -b is local the orbital time. For any location in 
the circumplanetary disk grains whose size obeys the latter re- 
lation are the fastest to migrate towards pressure maxima (i.e. 
towards the planet in this case given that pressure decreases 
monotonically from the planet outwards). Much smaller par- 
ticles will be entrained with the gas and migrate slowlier to the 
planet, while much larger boulders will follow perturbed kep- 
lerian orbits around the planet (Weidenschilling 1977). For 
isothermal and FLD runs the fastest migrating grains have 
sizes in the range 0.2 to 0.8 m, and the stopping time is al- 
most the same in both cases. We note that our calculation is 
strictly valid only down to 0. 15 AU from the planet; at smaller 
distances the assumption of the Epstein regime breaks down. 
For the adiabatic case, the Epstein regime is valid for bodies 
up to 10 m in size because the density of the circumplanetary 
envelope is very low. However, since the gas there is very hot, 
the fastest migrating grains are smaller than for the isothermal 
and FLD cases, having sizes in the range 0.06 to 0.1 m. 

Based on these results, we can expect that only very large 
boulders, i.e. those with sizes comfortably larger than indi- 
cated by the curves in Figure 4, e.g. in the range > 10- 100 m, 
may experience a slow enough radial decay to allow the as- 
sembly of satellites in the disk. However, addressing whether 



or not satellites can really form will require to include grain 
growth and the Stokes drag regime, as well as self-consistent 
initial sizes and velocities for boulders entering the circum- 
planetary disk. 



4. CONCLUSIONS 

We have shown the disk thermodynamics has an impact on 
all the most important aspects of disk-planet interaction for 
Jupiter-mass planets. The way migration and accretion are 
affected is qualitatively consistent with the findings by PM06 
for low-mass planets, although in that work the effect was 
quantitatively much stronger. 

Migration is slowed down relative to isothermal runs be- 
cause heating modifies the relative strength of the inner and 
outer spiral arm, or stifles the spiral perturbations completely 
as in the case of adiabatic runs. This effect partially counter- 
balances the faster migration found in our isothermal disks 
relative to standard estimates that neglect viscosity, self- 
gravity and the motion of the planet. If the slower migration 
seen in radiative transfer runs persists over > 100 orbits one 
would expect the mean migration rate to differ from the esti- 
mates of standard type II migration by no more than a factor of 
2, confirming the earlier conclusions reached by d'Angelo et 
al. (2003) using 2D simulations. The mass transport towards 
the planet is hampered by the increasing pressure gradient in 
non-isothermal runs. This will have to be investigated fur- 
ther by including a proper accretion model in the simulations 
(d'Angelo et al. 2003; KK06). 

The different structure of the circumplanetary envelope de- 
pending on disk thermodynamics that we find is qualitatively 
in agreement with the results of KK06, although our disks 
with flux-limited diffusion are less sub-keplerian than theirs. 
Satellite formation is not easy in the thick circumplanetary 
disk, but it may still be possible provided that the solid com- 
ponent is dominated by large boulders. 
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